-
Notifications
You must be signed in to change notification settings - Fork 238
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Refactor the Tian Approximation in the Reactive Fluid Transport Model #6161
base: main
Are you sure you want to change the base?
Refactor the Tian Approximation in the Reactive Fluid Transport Model #6161
Conversation
af3a3e0
to
f10d604
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sorry for the delayed review. This looks generally good, I just have a few comments. Please add a changelog entry that mentions this plugin is now easier to include in material models.
#ifndef _aspect_material_reaction_melt_tian2019_solubility_h | ||
#define _aspect_material_reaction_melt_tian2019_solubility_h |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I know you took this from the other reaction model, but that one sets a bad example here. These lines should correspond as closely as possible to the folder structure and file name. This is not technically necessary, but a convention we follow to make them easy to remember:
#ifndef _aspect_material_reaction_melt_tian2019_solubility_h | |
#define _aspect_material_reaction_melt_tian2019_solubility_h | |
#ifndef _aspect_material_model_reaction_model_tian2019_solubility_h | |
#define _aspect_material_model_reaction_model_tian2019_solubility_h |
|
||
template <int dim> | ||
Tian2019Solubility<dim>::Tian2019Solubility() | ||
= default; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If you make the constructor the default constructor then you usually do not need to implement it at all. Just remove it from here and from the header file.
/** | ||
* Declare the parameters this function takes through input files. | ||
*/ | ||
static | ||
void | ||
declare_parameters (ParameterHandler &prm); | ||
|
||
/** | ||
* Read the parameters from the parameter file. | ||
*/ | ||
void | ||
parse_parameters (ParameterHandler &prm); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Please put these functions at the end of the list of public functions here and in the source file. This is just a convention we (try to) follow, but it makes navigating the source files easier.
A typical class structure in ASPECT looks like this (with many functions being optional):
class ClassName
{
constructor();
initialize();
all_other_minor_functions();
main_function(in this case melt_fraction);
declare_parameters();
parse_parameters();
private:
... all member variables and private functions
}
|
||
|
||
/** | ||
* Compute the free fluid fraction that can be present in the material based on the |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
can be present
or is present
? I think you mean is present
because can be present
would be 100%?
/** | ||
* Compute the free fluid fraction that can be present in the material based on the | ||
* fluid content of the material and the fluid solubility for the given input conditions. | ||
* @p in and @p melt_fractions need to have the same size. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
there is no melt_fractions
variable
* @param q unsigned int from 0-3 indexing which rock phase the equilbrium | ||
* bound water content is being calculated for |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Oh, is that really the meaning of q
? We usually use q
as an index for the quadrature point to index into the vectors that are contained inside in
. Please clarify and if you use this variable for something else, maybe find a different name unless it has to be q
for other reasons.
/** | ||
* Parameters for the solubility model of Tian et al., 2019. | ||
*/ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think this comment here is redundant now. All of these parameters are from Tian et al since this is the Tian et al plugin.
std::vector<double> LR_peridotite_poly_coeffs {-19.0609, 168.983, -630.032, 1281.84, -1543.14, 1111.88, -459.142, 95.4143, 1.97246}; | ||
std::vector<double> csat_peridotite_poly_coeffs {0.00115628, 2.42179}; | ||
std::vector<double> Td_peridotite_poly_coeffs {-15.4627, 94.9716, 636.603}; | ||
|
||
std::vector<double> LR_gabbro_poly_coeffs {-1.81745, 7.67198, -10.8507, 5.09329, 8.14519}; | ||
std::vector<double> csat_gabbro_poly_coeffs {-0.0176673, 0.0893044, 1.52732}; | ||
std::vector<double> Td_gabbro_poly_coeffs {-1.72277, 20.5898, 637.517}; | ||
|
||
std::vector<double> LR_MORB_poly_coeffs {-1.78177, 7.50871, -10.4840, 5.19725, 7.96365}; | ||
std::vector<double> csat_MORB_poly_coeffs {0.0102725, -0.115390, 0.324452, 1.41588}; | ||
std::vector<double> Td_MORB_poly_coeffs {-3.81280, 22.7809, 638.049}; | ||
|
||
std::vector<double> LR_sediment_poly_coeffs {-2.03283, 10.8186, -21.2119, 18.3351, -6.48711, 8.32459}; | ||
std::vector<double> csat_sediment_poly_coeffs {-0.150662, 0.301807, 1.01867}; | ||
std::vector<double> Td_sediment_poly_coeffs {2.83277, -24.7593, 85.9090, 524.898}; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
All of these numbers are fixed and will never change during the program runtime. You can tell this to the compiler, which will then optimize the code further and also prevent you from accidentally changing the values by changing their type from std::vector<double>
to constexpr std::array<double,N>
where N is the number of values in this array. Initialization and access to these arrays can remain unchanged.
* and 50 GPa for sediment. These cutoff pressures were determined by extending the pressure range in Tian et al. (2019) | ||
* and observing where the maximum allowed water contents jump towards infinite values. | ||
*/ | ||
const std::vector<double> pressure_cutoffs {10, 26, 16, 50}; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
same here, make this a constexpr std::array
prm.enter_subsection("Tian 2019 model"); | ||
{ | ||
tian2019_model.initialize_simulator (this->get_simulator()); | ||
tian2019_model.parse_parameters(prm); | ||
} | ||
prm.leave_subsection(); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
So you put the parameters of the Tian model in their own subsection now. However, this means the parameter files of everyone using the reactive fluid flow model with the 3.0 version will need to be updated. Please include a changelog entry that mentions this. Additionally, you could extend our update scripts in contrib/utilities/update_scripts/prm_files/
to automatically move this parameter. The closest operation that does something like this is probably this function. However, the update process looks a bit cryptic, I leave this up to you to decide if this is important enough.
I was looking through the Reactive fluid transport model code and saw how the katz 2003 model is in it's own reaction model. I think I recall at the hackathon that the ultimate goal for this material model was to expand it by creating new reaction models, so I went and refactored the tian approximation into it's own reaction.